Dataset Exploratory Analysis

Module-Trait associations

datTraits = datMeta %>% dplyr::select(Diagnosis_, Region, Sex, Age, PMI, RNAExtractionBatch) %>%
            rename('Diagnosis' = Diagnosis_, 'ExtractionBatch' = RNAExtractionBatch)

# Recalculate MEs with color labels
ME_object = datExpr %>% t %>% moduleEigengenes(colors = dataset$Module)
MEs = orderMEs(ME_object$eigengenes)

# Calculate correlation between eigengenes and the traits and their p-values
moduleTraitCor = MEs %>% apply(2, function(x) hetcor(x, datTraits, ML=TRUE)$correlations[1,-1]) %>% t
rownames(moduleTraitCor) = colnames(MEs)
colnames(moduleTraitCor) = colnames(datTraits)
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nrow(datExpr))

# In case there are any NAs
if(sum(!complete.cases(moduleTraitCor))>0){
  print(paste0(sum(is.na(moduleTraitCor)),' correlation(s) could not be calculated')) 
}

## [1] "1 correlation(s) could not be calculated"

moduleTraitCor = moduleTraitCor[complete.cases(moduleTraitCor),]
moduleTraitPvalue = moduleTraitPvalue[complete.cases(moduleTraitCor),]

# Sort moduleTraitCor by Diagnosis
moduleTraitCor = moduleTraitCor[order(moduleTraitCor[,1], decreasing=TRUE),]
moduleTraitPvalue = moduleTraitPvalue[order(moduleTraitCor[,1], decreasing=TRUE),]

# Create text matrix for the Heatmap
textMatrix = paste0(signif(moduleTraitCor, 2), ' (', signif(moduleTraitPvalue, 1), ')')
dim(textMatrix) = dim(moduleTraitCor)


labeledHeatmap(Matrix = moduleTraitCor, xLabels = names(datTraits), yLabels =  gsub('ME','',rownames(moduleTraitCor)), 
               yColorWidth=0, colors = brewer.pal(11,'PiYG'), bg.lab.y = gsub('ME','',rownames(moduleTraitCor)),
               textMatrix = textMatrix, setStdMargins = FALSE, cex.text = 0.8, cex.lab.y = 0.75, zlim = c(-1,1),
               main = paste('Module-Trait relationships'))

rm(ME_object, datTraits, MEs, moduleTraitCor, moduleTraitPvalue, textMatrix)

Selecting the three modules with highest (absolute) correlation to Diagnosis: #C17FFF #00C094 and #C49A00

pca = datExpr %>% prcomp

plot_data = data.frame('ID'=rownames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2]) %>%
            left_join(dataset, by='ID') %>% dplyr::select(ID, PC1, PC2, Module, gene.score) %>%
            mutate(ImportantModules = ifelse(Module=='#C17FFF', '#C17FFF',
                                      ifelse(Module=='#00C094', '#00C094',
                                      ifelse(Module=='#C49A00', '#C49A00', 'Others')))) %>%
            mutate(color = ifelse(ImportantModules=='Others','gray',ImportantModules),
                   alpha = ifelse(ImportantModules=='Others', 0.1, 0.5))

ggplotly(plot_data %>% ggplot(aes(PC1, PC2, color=ImportantModules)) + 
         geom_point(alpha=plot_data$alpha, color=plot_data$color) + theme_minimal() + 
           ggtitle('Modules with strongest relation to Diagnosis'))

SFARI genes and Module-Diagnosis significance

Modules with the strongest module-diagnosis correlation should have the highest percentage of SFARI Genes, but the plot shows this relation in modules with correlations close to 0, and thi behaviour is the same for the other clusterings I have tried…

plot_data = dataset %>% mutate('hasSFARIscore' = gene.score!='None') %>% 
            group_by(Module, MTcor, hasSFARIscore) %>% summarise(p=n()) %>% 
            left_join(dataset %>% group_by(Module) %>% summarise(n=n()), by='Module') %>% 
            mutate(p=round(p/n*100,2)) 

for(i in 1:nrow(plot_data)){
  this_row = plot_data[i,]
  if(this_row$hasSFARIscore==FALSE & this_row$p==100){
    new_row = this_row
    new_row$hasSFARIscore = TRUE
    new_row$p = 0
    plot_data = plot_data %>% rbind(new_row)
  }
}

plot_data = plot_data %>% filter(hasSFARIscore==TRUE)

ggplotly(plot_data %>% ggplot(aes(MTcor, p, size=n)) + geom_smooth(color='gray', se=FALSE) +
         geom_point(color=plot_data$Module, aes(id=Module)) +
         xlab('Module-Diagnosis correlation') + ylab('% of SFARI genes') +
         theme_minimal() + theme(legend.position = 'none'))

rm(i, this_row, new_row)

Gene Significance

Gene significance close to 0 means the gene is not significant and negative values means the gene is underexpressed in Autism samples. So the important thing about gene significance is it’s absolute value

plot_data = data.frame('ID'=rownames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2]) %>%
            left_join(dataset, by='ID')

ggplotly(plot_data %>% ggplot(aes(PC1, PC2, color=GS)) + geom_point(alpha=0.5) + 
         scale_color_gradientn(colours=c('#F8766D','white','#00BFC4')) +
  theme_minimal() + ggtitle('Gene Significance'))

Module Membership and Gene Significance relation

plot_data = dataset %>% dplyr::select(MM.C17FFF, GS, gene.score) %>% filter(dataset$Module=='#C17FFF')

ggplotly(plot_data %>% ggplot(aes(MM.C17FFF, GS, color=gene.score)) + geom_point(alpha=0.5) + 
         scale_color_manual(values=SFARI_colour_hue(r=c(1:6,8,7))) + 
         theme_minimal())

plot_data = dataset %>% dplyr::select(MM.C49A00, GS, gene.score) %>% filter(dataset$Module=='#C49A00')

ggplotly(plot_data %>% ggplot(aes(MM.C49A00, GS, color=gene.score)) + geom_point(alpha=0.5) + 
         scale_color_manual(values=SFARI_colour_hue(r=c(2:6,8,7))) + 
         theme_minimal())

plot_data = dataset %>% dplyr::select(MM.00C094, GS, gene.score) %>% filter(dataset$Module=='#00C094')


ggplotly(plot_data %>% ggplot(aes(MM.00C094, GS, color=gene.score)) + geom_point(alpha=0.5) + 
         scale_color_manual(values=SFARI_colour_hue(r=c(3:6,8,7))) + 
         theme_minimal())

Identifying important genes

Selecting the modules with the highest correlation to Diagnosis, and, from them, the genes with the highest module membership-(absolute) gene significance

*Ordered by \(\frac{MM+GS}{2}\)

Only one SFARI Gene in the two lists …

top_genes_1 = dataset %>% dplyr::select(ID, MM.C17FFF, GS, gene.score) %>% filter(dataset$Module=='#C17FFF') %>%
              rename('MM' = 'MM.C17FFF') %>% mutate(importance = (MM+abs(GS))/2) %>% arrange(by=-importance) %>%
              top_n(10)

kable(top_genes_1, caption='Top 10 genes for module #C17FFF')

top_genes_2 = dataset %>% dplyr::select(ID, MM.C49A00, GS, gene.score) %>% filter(dataset$Module=='#C49A00') %>%
              rename('MM' = 'MM.C49A00') %>% mutate(importance = (MM+abs(GS))/2) %>% arrange(by=-importance) %>%
              top_n(10)

kable(top_genes_2, caption='Top 10 genes for module #C49A00')

top_genes_3 = dataset %>% dplyr::select(ID, MM.00C094, GS, gene.score) %>% filter(dataset$Module=='#00C094') %>%
              rename('MM' = 'MM.00C094') %>% mutate(importance = (MM+abs(GS))/2) %>% arrange(by=-importance) %>%
              top_n(10)

kable(top_genes_3, caption='Top 10 genes for module #00C094')

plot_data = data.frame('ID'=rownames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2]) %>%
            left_join(dataset, by='ID') %>% dplyr::select(ID, PC1, PC2, Module, gene.score) %>%
            mutate(color = ifelse(Module=='#C17FFF', '#C17FFF',
                           ifelse(Module=='#00C094', '#00C094',
                           ifelse(Module=='#C49A00', '#C49A00', 'gray')))) %>%
            mutate(alpha = ifelse(color %in% c('#C17FFF', '#00C094', '#C49A00') & 
                                  ID %in% c(as.character(top_genes_1$ID), 
                                            as.character(top_genes_2$ID),
                                            as.character(top_genes_3$ID)), 1, 0.1))

plot_data %>% ggplot(aes(PC1, PC2)) + geom_point(alpha=plot_data$alpha, color=plot_data$color) + 
              theme_minimal() + ggtitle('Important genes identified through WGCNA')